home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / src_original / zhpmv.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  8.5 KB  |  274 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE ZHPMV ( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY )
  5. *     .. Scalar Arguments ..
  6.       COMPLEX*16         ALPHA, BETA
  7.       INTEGER            INCX, INCY, N
  8.       CHARACTER*1        UPLO
  9. *     .. Array Arguments ..
  10.       COMPLEX*16         AP( * ), X( * ), Y( * )
  11. *     ..
  12. *
  13. *  Purpose
  14. *  =======
  15. *
  16. *  ZHPMV  performs the matrix-vector operation
  17. *
  18. *     y := alpha*A*x + beta*y,
  19. *
  20. *  where alpha and beta are scalars, x and y are n element vectors and
  21. *  A is an n by n hermitian matrix, supplied in packed form.
  22. *
  23. *  Parameters
  24. *  ==========
  25. *
  26. *  UPLO   - CHARACTER*1.
  27. *           On entry, UPLO specifies whether the upper or lower
  28. *           triangular part of the matrix A is supplied in the packed
  29. *           array AP as follows:
  30. *
  31. *              UPLO = 'U' or 'u'   The upper triangular part of A is
  32. *                                  supplied in AP.
  33. *
  34. *              UPLO = 'L' or 'l'   The lower triangular part of A is
  35. *                                  supplied in AP.
  36. *
  37. *           Unchanged on exit.
  38. *
  39. *  N      - INTEGER.
  40. *           On entry, N specifies the order of the matrix A.
  41. *           N must be at least zero.
  42. *           Unchanged on exit.
  43. *
  44. *  ALPHA  - COMPLEX*16      .
  45. *           On entry, ALPHA specifies the scalar alpha.
  46. *           Unchanged on exit.
  47. *
  48. *  AP     - COMPLEX*16       array of DIMENSION at least
  49. *           ( ( n*( n + 1 ) )/2 ).
  50. *           Before entry with UPLO = 'U' or 'u', the array AP must
  51. *           contain the upper triangular part of the hermitian matrix
  52. *           packed sequentially, column by column, so that AP( 1 )
  53. *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
  54. *           and a( 2, 2 ) respectively, and so on.
  55. *           Before entry with UPLO = 'L' or 'l', the array AP must
  56. *           contain the lower triangular part of the hermitian matrix
  57. *           packed sequentially, column by column, so that AP( 1 )
  58. *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
  59. *           and a( 3, 1 ) respectively, and so on.
  60. *           Note that the imaginary parts of the diagonal elements need
  61. *           not be set and are assumed to be zero.
  62. *           Unchanged on exit.
  63. *
  64. *  X      - COMPLEX*16       array of dimension at least
  65. *           ( 1 + ( n - 1 )*abs( INCX ) ).
  66. *           Before entry, the incremented array X must contain the n
  67. *           element vector x.
  68. *           Unchanged on exit.
  69. *
  70. *  INCX   - INTEGER.
  71. *           On entry, INCX specifies the increment for the elements of
  72. *           X. INCX must not be zero.
  73. *           Unchanged on exit.
  74. *
  75. *  BETA   - COMPLEX*16      .
  76. *           On entry, BETA specifies the scalar beta. When BETA is
  77. *           supplied as zero then Y need not be set on input.
  78. *           Unchanged on exit.
  79. *
  80. *  Y      - COMPLEX*16       array of dimension at least
  81. *           ( 1 + ( n - 1 )*abs( INCY ) ).
  82. *           Before entry, the incremented array Y must contain the n
  83. *           element vector y. On exit, Y is overwritten by the updated
  84. *           vector y.
  85. *
  86. *  INCY   - INTEGER.
  87. *           On entry, INCY specifies the increment for the elements of
  88. *           Y. INCY must not be zero.
  89. *           Unchanged on exit.
  90. *
  91. *
  92. *  Level 2 Blas routine.
  93. *
  94. *  -- Written on 22-October-1986.
  95. *     Jack Dongarra, Argonne National Lab.
  96. *     Jeremy Du Croz, Nag Central Office.
  97. *     Sven Hammarling, Nag Central Office.
  98. *     Richard Hanson, Sandia National Labs.
  99. *
  100. *
  101. *     .. Parameters ..
  102.       COMPLEX*16         ONE
  103.       PARAMETER        ( ONE  = ( 1.0D+0, 0.0D+0 ) )
  104.       COMPLEX*16         ZERO
  105.       PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
  106. *     .. Local Scalars ..
  107.       COMPLEX*16         TEMP1, TEMP2
  108.       INTEGER            I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY
  109. *     .. External Functions ..
  110.       LOGICAL            LSAME
  111.       EXTERNAL           LSAME
  112. *     .. External Subroutines ..
  113.       EXTERNAL           XERBLA
  114. *     .. Intrinsic Functions ..
  115.       INTRINSIC          DCONJG, DBLE
  116. *     ..
  117. *     .. Executable Statements ..
  118. *
  119. *     Test the input parameters.
  120. *
  121.       INFO = 0
  122.       IF     ( .NOT.LSAME( UPLO, 'U' ).AND.
  123.      $         .NOT.LSAME( UPLO, 'L' )      )THEN
  124.          INFO = 1
  125.       ELSE IF( N.LT.0 )THEN
  126.          INFO = 2
  127.       ELSE IF( INCX.EQ.0 )THEN
  128.          INFO = 6
  129.       ELSE IF( INCY.EQ.0 )THEN
  130.          INFO = 9
  131.       END IF
  132.       IF( INFO.NE.0 )THEN
  133.          CALL XERBLA( 'ZHPMV ', INFO )
  134.          RETURN
  135.       END IF
  136. *
  137. *     Quick return if possible.
  138. *
  139.       IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
  140.      $   RETURN
  141. *
  142. *     Set up the start points in  X  and  Y.
  143. *
  144.       IF( INCX.GT.0 )THEN
  145.          KX = 1
  146.       ELSE
  147.          KX = 1 - ( N - 1 )*INCX
  148.       END IF
  149.       IF( INCY.GT.0 )THEN
  150.          KY = 1
  151.       ELSE
  152.          KY = 1 - ( N - 1 )*INCY
  153.       END IF
  154. *
  155. *     Start the operations. In this version the elements of the array AP
  156. *     are accessed sequentially with one pass through AP.
  157. *
  158. *     First form  y := beta*y.
  159. *
  160.       IF( BETA.NE.ONE )THEN
  161.          IF( INCY.EQ.1 )THEN
  162.             IF( BETA.EQ.ZERO )THEN
  163.                DO 10, I = 1, N
  164.                   Y( I ) = ZERO
  165.    10          CONTINUE
  166.             ELSE
  167.                DO 20, I = 1, N
  168.                   Y( I ) = BETA*Y( I )
  169.    20          CONTINUE
  170.             END IF
  171.          ELSE
  172.             IY = KY
  173.             IF( BETA.EQ.ZERO )THEN
  174.                DO 30, I = 1, N
  175.                   Y( IY ) = ZERO
  176.                   IY      = IY   + INCY
  177.    30          CONTINUE
  178.             ELSE
  179.                DO 40, I = 1, N
  180.                   Y( IY ) = BETA*Y( IY )
  181.                   IY      = IY           + INCY
  182.    40          CONTINUE
  183.             END IF
  184.          END IF
  185.       END IF
  186.       IF( ALPHA.EQ.ZERO )
  187.      $   RETURN
  188.       KK = 1
  189.       IF( LSAME( UPLO, 'U' ) )THEN
  190. *
  191. *        Form  y  when AP contains the upper triangle.
  192. *
  193.          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  194.             DO 60, J = 1, N
  195.                TEMP1 = ALPHA*X( J )
  196.                TEMP2 = ZERO
  197.                K     = KK
  198.                DO 50, I = 1, J - 1
  199.                   Y( I ) = Y( I ) + TEMP1*AP( K )
  200.                   TEMP2  = TEMP2  + DCONJG( AP( K ) )*X( I )
  201.                   K      = K      + 1
  202.    50          CONTINUE
  203.                Y( J ) = Y( J ) + TEMP1*DBLE( AP( KK + J - 1 ) )
  204.      $                         + ALPHA*TEMP2
  205.                KK     = KK     + J
  206.    60       CONTINUE
  207.          ELSE
  208.             JX = KX
  209.             JY = KY
  210.             DO 80, J = 1, N
  211.                TEMP1 = ALPHA*X( JX )
  212.                TEMP2 = ZERO
  213.                IX    = KX
  214.                IY    = KY
  215.                DO 70, K = KK, KK + J - 2
  216.                   Y( IY ) = Y( IY ) + TEMP1*AP( K )
  217.                   TEMP2   = TEMP2   + DCONJG( AP( K ) )*X( IX )
  218.                   IX      = IX      + INCX
  219.                   IY      = IY      + INCY
  220.    70          CONTINUE
  221.                Y( JY ) = Y( JY ) + TEMP1*DBLE( AP( KK + J - 1 ) )
  222.      $                           + ALPHA*TEMP2
  223.                JX      = JX      + INCX
  224.                JY      = JY      + INCY
  225.                KK      = KK      + J
  226.    80       CONTINUE
  227.          END IF
  228.       ELSE
  229. *
  230. *        Form  y  when AP contains the lower triangle.
  231. *
  232.          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  233.             DO 100, J = 1, N
  234.                TEMP1  = ALPHA*X( J )
  235.                TEMP2  = ZERO
  236.                Y( J ) = Y( J ) + TEMP1*DBLE( AP( KK ) )
  237.                K      = KK     + 1
  238.                DO 90, I = J + 1, N
  239.                   Y( I ) = Y( I ) + TEMP1*AP( K )
  240.                   TEMP2  = TEMP2  + DCONJG( AP( K ) )*X( I )
  241.                   K      = K      + 1
  242.    90          CONTINUE
  243.                Y( J ) = Y( J ) + ALPHA*TEMP2
  244.                KK     = KK     + ( N - J + 1 )
  245.   100       CONTINUE
  246.          ELSE
  247.             JX = KX
  248.             JY = KY
  249.             DO 120, J = 1, N
  250.                TEMP1   = ALPHA*X( JX )
  251.                TEMP2   = ZERO
  252.                Y( JY ) = Y( JY ) + TEMP1*DBLE( AP( KK ) )
  253.                IX      = JX
  254.                IY      = JY
  255.                DO 110, K = KK + 1, KK + N - J
  256.                   IX      = IX      + INCX
  257.                   IY      = IY      + INCY
  258.                   Y( IY ) = Y( IY ) + TEMP1*AP( K )
  259.                   TEMP2   = TEMP2   + DCONJG( AP( K ) )*X( IX )
  260.   110          CONTINUE
  261.                Y( JY ) = Y( JY ) + ALPHA*TEMP2
  262.                JX      = JX      + INCX
  263.                JY      = JY      + INCY
  264.                KK      = KK      + ( N - J + 1 )
  265.   120       CONTINUE
  266.          END IF
  267.       END IF
  268. *
  269.       RETURN
  270. *
  271. *     End of ZHPMV .
  272. *
  273.       END
  274.